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1 Introduction 

1.1 The problem and background 

1.1.1 The problem 

Let (Q, T, P) be a probability space, B°(f) = t, and (B 1 (t), . . . , B d (t)) be a d-dimensional 
standard Brownian motion. C™(R N ; R N ) denotes the set of R N -valued infinitely dif- 
ferentiable functions defined in R N whose derivatives are all bounded. Our interest 
is in weak approximation, that is to say, approximation of (Pf/)(x) = E[f(X(t,x))] 
where / e C^(R N ;R) and X(t, x) is a solution to the stochastic differential equation 
written in the Stratonovich form: 

i m 

X(t,x) = x+Y I Vi(X(s,x)) o dB\s), (1.1) 

where V, e q°(R N ; R N ) for i = 0, 1, . . . , d. Here, V 7 ; e C£° (r n ; R n ) is considered to be 
a vector field in the following way: 

Vif(x) = Yj v lMdx~W> for / e c r( RN ; R )- 

7=1 ' 

It is well-known (e. g. 1 13 1) that E[f (X(l, x))] is equal to u(l, x) where u is the solution 
to the following partial differential equation for L — V + (1/2) Vf : 

du 

—(t,x) = Lu, w(0, x) = /(*). (1.2) 

1.1.2 Background 

A number of studies on numerical calculations of this problem have been conducted 
as there is a great demand for it in various fields. One often encounters this type 
of calculation particularly in mathematical finance. For example, the price of a 
financial derivative written on the diffusion X(t, x) is obtained by the calculation of 
E[f(X(T,x))]. 

There are two approaches to the problem: PDE approach and simulation. The 
former one involves solving the partial differential equation dl.2t numerically. This 
method works only when L is elliptic and the dimension is relatively small. We do 
not go into details on the subject here but refer to |19|. These conditions are not 
necessarily satisfied in many practical problems so we are forced to take the other 
approach which is called the probabilistic method or simulation. In this paper, we 
focus on this approach. 

Usually, the Euler-Maruyama scheme is used to discretize X(t, x) during simu- 
lations to weakly approximate X(t,x). It is shown in |18|, |23[, |24[, and |28[ that the 
new higher-order scheme introduced by Kusuoka in 1151 calculates some finance 
problems much faster than the Euler-Maruyama scheme. Lyons and Victoir exten- 
sively developed the scheme in |21| using the notion of free Lie algebra. Recent 
developments can be found in |2 [ and 1 10 [. 

We will discuss the reason why higher order schemes greatly improve the speed 
of numerical weak approximation in the later part of this paper (Section|6]l. 
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2.2.3 Our results 

In this paper, we describe how successfully we constructed in Theorem 11.31 and 
Corollary 1 1.41 a new higher order weak approximation scheme for a broad class of 
stochastic differential equations. This scheme owes a great deal to the scheme shown 
in [15;| and to the cubature method on Wiener space introduced in [21 J . 

An intuitive explanation of the scheme is as follows. We construct the ODE (or- 
dinary differential equation)-valued random variable whose average approximates 
the given stochastic differential equation. From this random variable, an ODE itself 
is able to be drawn at one time. 

This scheme has a remarkable advantage that once an ODE is drawn, the con- 
ventional Runge-Kutta method can be applied so as to approximate the ODE. The 
approximating random variable is constructed using Theorem ll.3l and Theorem ll.6l 
and can be approximated by the Runge-Kutta method for ODEs via Theorem l4.15l 

We should note that another higher-order weak approximation method is intro- 
duced in 1 25 1 . Although the algorithm in |25[ and the new method presented in this 
paper are based on the same scheme (|15| 121 J) and have many common features, 
algorithms themselves differ significantly. 



1.2 Notation 

Let A = {vq, V\, . . . , Vj) be an alphabet where d e Z>i and A* denote the set of all 
words consisting of the elements of A. The empty word 1 is the identity of A*. For u = 
■ • • Vj n e A*, \u\ and ||u|| are defined by \u\ = n and ||w|| = \u\ + card({fc| ^ = 0}) where 
card(S) denotes the cardinality of a set S. Here, || • || is related to the scaling property 
of the Brownian motion. A*„ and A* <m denote [w € A*\ \w\ = m\ and {w e A*\ \zv\ < m), 
respectively. Let R(A) be the R-coefficient free algebra with basis A* and R((A)) be 
the set of all R-coefficient formal series with basis A*. Then, R(A) is a sub-R-algebra 
of R((A)). We call an element of R(A) a non-commutative polynomial. P e R((A)) 
is written as 

P = ^ (P, w) w or "S\ a w w, 



weA' 



weA' 



where (P, w) = a w e R denotes the coefficient of w. Let 

R(A) m = {P e R<A) | (P r w) = 0, if ||w|| ± m) 
The algebra structure is defined as usual, i.e. 



\weA' 



weA" 



The Lie bracket is defined as [x, y] = xy— yx for x,y e R((A)). For w = ■ ■ ■ Vj n e A*, 
l(w) denotes [vj u [vj 2 , [. . . , [vi l: _ u z>,J . . . ]]]. We define Xk (A) as the set of Lie polyno- 
mials in R(A) and Xr((A)) as the set of Lie series. This means that -£r(A) is the 
smallest sub-R-module of R(A) including A and is closed under the Lie bracket, 
and that Xr((A)) is the set of elements of R((A)) whose homogeneous components 
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belong to -£r(A). We note that Lie polynomials correspond to vector fields while 
general polynomials do not necessarily. For m e Z> , let j m be a map defined by 

a w w = ^ a w w. 

\weA' ) \\w\\<m 

For arbitrary P, Q e R(A), the inner product (P, Q) is defined by 

(P,Q)= J^frwXQM- 

zveA" 

Moreover we let ||P|| 2 = «P,P» 1/2 for P e R<A>. For P e R«A» with (P,l) = 
0, we can define exp(P) as 1 + Y*kL\P k /kl. In addition, log(Q) can be defined as 
EfcLi(-l)* _1 (Q ~ lf/k for Q £ R«A» with (Q, 1) = 1. Then the following relations 
hold: 

log(exp(P)) = P and exp(log(Q)) = Q. 

By the natural identification R«A» « R°°, we can induce the direct product 
topology into R((A)). Then, R((A)) becomes a Polish space by the topology. We can 
also consider its Borel a-algebra S(R((A))), R((A))-valued random variables, their 
expectations, and other notions as usual. 

Let O be the homomorphism between R(A) and the R-algebra consisting of 
smooth differential operators over R N such that 

0(1) = Id, 

(1-3) 

®{v h ---Vi n ) = V k ---V in for /J,..., z„ e {0,1,. . 

Considering the scaling property of the Brownian motion, we define the rescaling 
operator W s depending on || • ||. For s e R >0 , W s : R«A» — > R«A» is defined as 
follows: 

w v 



Yj p ™ =Yj S "' /Zp ™ where P m eR<A)„ 



\m=0 ) m=0 

For a smooth vector field V, i. e. an element of CjJ° (R N ; R n ), exp (V) (x) denotes 
the solution at time 1 of the ordinary differential equation 

^ = V (z ( ) , z = x. 
We also define ||V||c« for V e q°(R N ;R N ) as follows: 
||V||=sup{|V(x)|; «R N ) 

= sup fy$(Uu U 2 ,..., U„)\ ;xeU N and |U,| = 1, for i = 1, . . . , n} 

IMIc"=£|M|- 

2=0 

Here denotes the fcth order total differential of V, i.e. 

N N N , 

vg ^ u») = x Z • • • E ^- .4. WLZ i' • • • u " ei 

i=1 /i=l ;„=1 ^ 

where each e, denotes an N-dimensional unit vector, {e\,..., e N ] forms an orthonor- 
mal basis of R N , and U ! k is the ;'th component of e R N . 



5 



1.3 Main results 

Since in this paper we deal with the operators that are not necessarily linear with 
respect to time t, we introduce the following definition: 

Definition 1.1 A map gfrom C~(]R N ;]R N ) to the set of all maps from R N to R N z's called 
an integration scheme of order m if there exists a positive constant C m such that 



sup \g(W)(x) - exp (W)(x)\ < C m \\W\\ 

xeR N 



m+1 

Qn+1 



(1.4) 



for all W G C"(R N ; R N ). Let IS(m) be the set of all integration schemes of order m. 

This definition is a generalization of the usual order of approximation. 

Definition 1.2 For Z\,z% e Xr((A)), we define Z2HZ1 as log(exp(z2)exp(zi)). Then from 
the definition, for Zi,z 2 ,z 3 e X.r{{A)), 

(ziHz 2 )hz 3 = log(exp(zi)exp(z 2 )exp(z 3 )) = ZiH (z 2 hz 3 ), 

and so we can write for Z\, . . . , z„ € Xr((A)) 

ZiHz 2 h • • • hz„ = log(exp(z!) • • • exp(z„)) . (1.5) 

We notice that z 2 HZi € Xr((A)) if z lr z 2 e Xr((A)) from the Baker-Campbell- 
Hausdorff formula(|4 [). 

The following are the main results. 

Theorem 1.3 Let m > 1, M > 2, and . .,Zm be X.u{{A))-valued random variables. 
Assume that Z\, ... , Zm satisfy the fallowings: 



Z, = j m Z{ fori- 1, . . . , M, 
E[||/ w Zi[|J<oo fori = l,...,M, 



exp 



c „,+i 



< 00 for any a > 0. 



(1.6) 
(1.7) 

(1.8) 



Then for p G [1, 00) and arbitrary gx,--.,gu e IS{m), there exists a positive constant C Mj m 
such that 



sup |# (<P (W s (Zi))) o • • • o g M (O (W s (Z M ))) (x) 

ieR N 

- exp (O (W s O'm (Z m h • ■ ■ HZi)))) (x) 



— C#lMS 



(m+l)/2 



(1.9) 



/or s e (0, 1] rofcere C m/ M depends only on m and M. Here for functions f and g, f o g(x) 
denotes f (g{x)) as usual. 
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For i = 1, . . . , d, and j — 1, . . . , M, let S'. be R-valued Gaussian random variables 
and for /, f = 1, . . . ,M, let e; and R;;< be real numbers such that 

M 

Y Cj = 1, E [S^] = 0, and E [sj-S^] = %6 !? (1.10) 
;'=i 

for i, i' = 1, . . . , d. We let S° = c ; for convenience. Taking (1.2) into account, we let 

Zi, . . . , Zm be random variables such that Zj - CjVg + Y,i=i $ % - v i f° r / = 1, . . . ,M and 
that 



E [j m (exp (Z : ) ■ ■ • exp (Z M ))] = /„ 



1 


1 


d \ 




exp 


^0 + 2 










;=i J 


) 



(1.11) 



In usual ODE cases, this type of approximation technique is known as a split- 
ting method (|12|). The stochastic versions of this technique are considered in |21 1 
and (251. 

Corollary 1.4 Suppose that the following UFG condition is satisfied: 
(UFG) There exist an integer I and (p U/U > e C£°(R N ; R) which satisfy 

<Z>(r(«))= Yj fu,wO(x(u')) (1.12) 

for any u e A* \ |l,i?o}- 
For j = l,...,MletZjbe X.u{(A))-valued random variables constructed as above and define 
linear operators Qi s \ for s e (0, 1] by 

(Q W /) (x) = E [/ (g (<P (W s (Zj))) o • • • o g((P (f s (Z M ))) (x))] (1.13) 

where f e q°(R N ; R) and g e IS{m). Then 

\\P s f - Q (S) /|L * Cs^ HgradCOL (1.14) 
w/iere C fs a positive constant. 

Remark 1.5 In jflTF , zf is shown that for the operator Q( s ) defined above, there exists a 
constant C and 

(P 5 f) (x) - (Q (s) /) (x) = Cs^ 1 + 0(s ( '" +3 > /2 ) 
holds. This means that the Romberg extrapolation can be applied to our new algorithm. 

The intuitive understanding is that once we find the random variables Zi, . . . , Zm, 
we can numerically approximate exp(Z,(a>)) by applying the integration scheme g, 
repeatedly for each i as seen in jl.9) in Theorem |1.3l Therefore, our primary interest 
is in finding Z\, . . . , Zm- 

Theorem 1.6 Let m — 5 and M — 2. Then dl.lH holds if and only if 

+ V2(2w-1) n V2(2"-l) D 
ci = , c 2 = 1 ± -i , R u = m 

, y/2 (2u - 1) 

R22 = 1 + M ± a/2 (2m - 1), R12 = -« + — — 

for some u > 1/2. 



(1.15) 
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Remark 1.7 We can show that in the case where m = 7 and M — 3 there is no solution 
to Pi) . 

Now that we have obtained the random variables satisfying dl.llb , we need a 
practical way of approximating these integration schemes g\, . . . , gM- We success- 
fully extend applicability of the general Runge-Kutta method to ODEs to find that 
it belongs to IS(m). 

Let A = (fl, ; ). ._ k with a,j e IR and b = f (£>i, . . . , b K ) e R K . If (A, b) satisfies the 

m-th-order conditions defined as 14.8b in Section|H the K-stage Runge-Kutta method 
of order m in the sense of |6| can be written as follows: 



Yi (W s) = y + s £ at) W (y 7 (W s)) , 

K 

Y(y ; W, s) = y + s£ b t W (Y,( W, s)) 



' 1 (1.16) 

K 



for W e q° (r n ;R n ), s e R >0 , and y e R N - Let g(W)(y ) be Y(y ; W !)■ We show 
that g belongs to IS{m) in Theorem l4.15l 

Remark 1.8 Our scheme is fundamentally different from the class of numerical methods 
sometimes referred to as stochastic Runge-Kutta methods (15]126][27J). 



2 Proof of Theorem H31 

We split the left-hand side of as 

|| sup \g! (O s (Zi)) o ■ ■ • o gM (0 S (Z M )) (x) 

- exp (® s (j m (Z m h • • • hZi))) (x) 



(2.1) 



< || sup |exp (O s (Z x )) o ••• o exp (® s (Z M )) (x) 

- exp (® s (j m (.Z M H--- HZ X ))) Wll^ 
+ || sup \ gl (O s (Zi)) o • • ■ o g M (O s (Z M )) (x) 

XER N 

- exp (<D S (Zj)) o ... o exp (O s (Z M )) (x)| ||„. 

Evaluation of each term of the right-hand side of 02.1b will be given by Lemma |Z6l 
or d2.14b in this section. 

Proposition 2.1 

(1) For any V e q°(R N ; R N ), / e C° 3 (R N ; R), x e R N and « > 1, 

/(exp(fV)(*)) = g ^ (V k f) (x) + ^ (V" +I f) (exp(sV)(x)) ds. (2.2) 



s 



(2) For all z e -£r((A)) and n,m>\, 

f (exp (0(j m z)) (x)) ((/„z)*) /) (x 

^r4wll°(^ z ) M+1 )/|l ■ ( Z3 ) 

(n + 1)! II v ' Woo 



sup 

xeR N 



/c=0 



Proof Since we have 



/ (exp(f V) (x)) = (Vf) (exp (f V) (x)) , 



from the Taylor expansion, we obtain l|2.2f l by integration by parts and l !2.3t can be 
derived from 12.2)1 . □ 



Lemma 2.2 For a/1 n>l, there exists a constant C„ > such that 

ll^(/» z )/[L ^ c » lb» z ll 2 ll^^collc— 

for all z e R«A» and f e C°°(R N ; R). 
Proof Let p m be a map such that 

CO 

M=0 



(2.4) 



where a is a multi-index, a„ e Cf^R^R), and D" = „ 1 °""' „ N . Then we have 

oX l ...oX N 

m 

<2>(a>) = J] p,.((P(H7)), 
!=i 

for e A* \ {1}. Since there exists a constant C ny > such that 
\\pimw))f\\ x <C w ,i sup HD^gradtf))^, 

ae(Z> ) N 
\a\=i-l 

we see that there exists a constant C', > such that 



l<IMI<n 



< E E c ^ K2 ' w)| su p ll Da (s rad (/))ll 

weA' i=l «e(Z> ) N 
1<|M|<« |a|=i-l 

<C„\\j n z\\ 2 sup llD^gradC/))^ 

ae(Z> ) N 
M<n-1 

^ C « l|/» 2 ll 2 l|g rad (/)|lc-i 

where C' n = card ({w e A* 1 1 < ||zo|| < n)) sup ^ (EK C w , t ). 

l<\\w\\<n 
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For simplicity of notation, we let <J> s (y) denote <t>(W s (y)) for an element y e 
-Cr((^)) m the following part. 

Lemma 2.3 (1) There exists a constant C m \ > such that 
sup |/ (exp (O s (; m z)) (*)) - (O s (/,„ exp (j m z)) f) (x)\ 



xeR N 



< C m , lS ^ 2 (l + ||7 m z|| 2 )"' +1 ||grad(/)|| c „ ( ,„ +1M (2.5) 



for zeXr((A)). 
(2) TTzere exisis a constant C M/ m > wfere M e Z>2 swc/i f/zaf 

sup |/ (exp (<D S (j m ((/ m z M ) H • • • H (;' m zi)))) (*)) 

xeR N 



- (O s (j m exp (;' m (0' m z M ) h • • ■ h (; m zi)))) /) (x)| 



< C m . M s< m+1 > /2 



1 + IMIz llgrad^llc^)-! ( 2 - 6 ) 



/orzi,...,z M 6 £r((A)). 
Proo/ From the fact that for z e X K ((A)) 

m 1 m 1 

/,„ (exp(/,„z)) = ^ - 0«Z)* - rf(/m(m+l) - im)((jmzf) 
k=0 ' *:=2 

and from 02.31 in Proposition ^. 11 we see that 
|/ (exp (<P(/ m z)) (*)) - (<P (/,„ (exp(; m z))) /) (x)| 

- f^TTv II* ( (; ™ 2)M+1 ) + E F (° " 7m) ( (7mZ) ")) ^ W 

V ; ' k=2 ' 

Since we have 

(JmZ) ~ ( /m(m+l) jm^j (/mZ) / 

the folio wings can be derived by applying Lemma [2.2l 

|/ (exp (CP (; mZ )) (x)) - (CP (;,„ (exp (; m z))) /) (*)| 



(2.7) 



m+l ^ 

- E F I* (0 m <'" +1 > ~ /») (^ mZ)fC ))4oo 
fc=2 ' °° 

m+l 

< C,„ ||(; m(m+ i) - /m)(/»z)*|| llgradWllc^D-i 



(2.8) 



C„,,i (l 



+ \\1mZ\ 



llgrad(/)|| 



C m(m+1)-1 



where C m and C m ,i are positive constants. Thus d2.5t is proved. 
Taking (/ m ZAt)n ■ ■ ■ H(/ m Zi) as z above and evaluating by 

m+l I m \ 

y, ||(/m(JB+l) ~~ /m) Om ((/m z M) H ■ • • H < C nl/ M 1 + ||/»wZ; 

fc=2 

we obtain l !2.6t . 



i=i 



m+l 
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Lemma 2.4 There exists a constant C M/ m > such that 



sup |/ (exp (<P S (/ M Zi)) o ■ • • o exp (<D S (/ m z M )) (x)) 



- (® s (jm exp (j m ((j m Z M ) H • • • H (/ m Zl)))) /) (x)\ 



M 

<C m ,Ms('» +1 )/ 2 £(l + ||y m z4)'" +1 ||grad(/)||, 



(2.9) 



1=1 



for z\, . . . ,Zm £ Xr((A)). Here C m ,M depends on m and M. 

Proof We prove the lemma by induction on M. When M = 1, j2.5b and J2.9t are 
equivalent. Assume that i2.9\ holds for M. Splitting the left-hand side of <2.9b for 
M + 1 as 

|/(exp (<D S (; m Zi)) o • • • o exp (<D S (; m z M +i)) (x)) 

- (<J> S (/,„ exp (; m ((;' m z M +i) H • • • H (; m zi)))) /) (x)| 
< |/ (exp (<D S (/ OT Zi)) o ■ • • o exp (<P S (/ m z M+ i)) (x)) 

- (<D S (j m exp (/ w ((/ m z M ) h • • ■ h (/ m zi)))) /) (exp (<P S (; m z M +i)) (x))| 
+ \(® s (jm exp (/„, ((7' m z M ) h • • • h (/ m zi)))) /) (exp (O s (j m z M +\)) (x)) 

- (<£ s (7m exp (;',„ ((/ m z M +i) H • • • H (; m zi)))) /) (x)| , 

we can apply the induction hypothesis and d2.5l l with & s (j m exp ((/,„Zm) h • • • h (j m z\))) f 
instead of / obtaining 

sup |/(exp (O s (/ m za)) o ■ • • o exp (<D S (/ m z M+1 )) (x)) 



(d> s (j,„ exp (j m ((/ m Z M +l) H • • • H (/mZl)))) /) (x)\ 




where C\ > is a constant depending on m and M. Hence, <2.9t is proved. 



□ 



From Lemma [2.3| and l2~4l we have the following result. 



Lemma 2.5 For all m > 1, ffere exzsfs a constant C M/ m > smc/z that 



sup |/ (exp (<P S 0' m zi)) o • • • o exp (<D S (j m z M )) (x)) 



- / (exp (<D S (;' m ((;' m z M )H • • • H(;' m Zi)))) (x)) 



M 



< C,„, M s (m+1)/2 £ (l + ||;'mz4) m+ ||grad(/) 



(2.10) 



i=i 



for all s e (0, 1], z a , . . . ,z M e £r((A)), and / e C°°(R N ; R). 
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Lemma 2.6 Let Z\,...,Zm be £.^((A))-valued random variables such that for m > 1, 
£ [||/mZj[| 2 J < oo for i = 1, . . . , M. Then, for p e [1, oo) there exists a constant C„ h u > 
such that 



sup |exp (cf> s (jmZi)) o ■ ■ • o exp (<P S (j m Z M )) (a 



xeR w 



■ exp (CD S (/« ({jmZ M )H ■ ■ ■ H(; m Z!)))) (x)\ 



< C,„, M s (m+1)/2 (2.11) 



/or any s e (0, 1]. 



Proof If for i e {1, ... , N\, f((x\ x N )) = x', then ] | grad^/) 1 1 C ,„ ( ,„ +M) _ a = 1 for all m > 1. 
Therefore, applying Lemma |Z51 for this /, we obtain 12.1H . □ 



(2.12) 



We note that in |29 1 a similar result to this Lemma is obtained. 

We now start discussion about the latter term of the right-hand side of (2.1). 

Proposition 2.7 There exists a constant C > such that 

|g(W)(x) - g(W)(y)\ < C ||W||* + i + |x - y| exp (||W|b) 
for g e IS{m) and W e q°(R N ; R N ). 
Proof Since Gronwall's inequality gives 

|exp (W) (x) - exp (W) (y)| < \x- y| exp (||W|| c i) , 
(2. 12) can be derived. 

Since g; e IS(m) and each Z; satisfies jl.8t , we see that for some Ci > 0, 

sup \g M (cP s (Z M )) W - exp (cp s (Z M )) (x) 

xeR n 

< ||C,„ ||cPs (2^)11^11^ < C lS ('» +1 )/ 2 . (2.13) 
From this fact and Proposition ^. 71 there exists a constant C4 > such that 



sup \g M -i (<£>s (Z M -i)) o g M s (Z M )) (x) - exp (<P S (Z M _i)) o exp (cp s (Z M )) W 

xeR N 



sup |g M -i (f s (Z M -i)) ° exp (<Ps (Z M )) W - exp (cp s (Z M -i)) o exp (cp s (Z M )) « 

xeR N 



sup |g M -i {®s (Z M -i)) ° gM {<I>s (Z M )) (x) - gM-i (®s (Z M -i)) ° exp (cp s (Z M )) (x) 

xeR N 



< C 2 s 



(m+l)/2 



C 3 ||cP s (Z M _i)||^i + sup \g u (CPs (Z M )) (x) - exp (cp s (Z M )) (x)\ exp (||(p s (Z M _ 



xeR N 



< C 4 s 



(m+l)/2 
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where Ci and C3 are positive constants. Inductively, 



sup \gi {® s (Zi)) o • ■ • o g M (0 S (Z M )) (x) 

xeR N 



exp (<D S (Zi)) o ■ ■ ■ o exp (<£ s (Z M )) (x) 



<C 5 s 



(m+l)/2 



(2.14) 



where C5 > 0. 

Lemma l2~6l and d2.14t complete the proof. 



3 Construction of the _CiR((A))-valued random variables Zi,..., Z M 
Lemma 3.1 For i= 1,...,M, let Y; be Gaussian random variables such that 
E [Yj] = and E [ Y, Y ; ] = R(i, j), for i, j = 1, . . . ,M 

where R(i, j) e R. Moreover, for i = 1, . . . ,M let nii e Z> be such that Yh=i m i zs even. 
Then we have 



e[y^--y-«] 



rll<i</<M (dij-J l<i<j<M 



\ d -i\l<i<i<M ee< - mi '" M > 

where e {m.\, ... , m M ) is a set of {dijh< ( -<f<M satisfying that du e Z>o and that 



Y d ji + 2d » + Y d 'i = m! 

l<;'<i i<j<M 



fori = l,... ,M. 



Proof Let Z = E^i We have 



e [ y T ,,,y m m ] = E 



(?Z mi • • • dz"' M 



exp 



( m ^ 
£>Y 



<9z mi • • • dz mM 



exp 



2 

\ \<i,j<M 



2'/ 2 (//2)! dz^-'-dz' 



z=0 
\Z/2 



^ R(i,j)ziZj 

J<i,j<M 



z=0 



where z = (zi,...,z M ) e R M . 
Let 



1< i</<M 



rff; G Z> and ^ d,y = - I . 

\<i<j<M J 



(3-2) 



(3.3) 
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Then 



\U2 



M 



\l/2 



£ R(i,j)ZiZj = + 2 £ R{i,])z iZj \ 

\<i,j<M ) \ !=1 \<i<j<M ) 



L 



{1/2V- 



M 



(d;,h<;<,<Mefi 1 il<i<;'<M ^"'/-j <=1 l<i</"<M 



E 



(//2)i 



(d,,h<;<,<Mefi ril<i<;<M (^';') 

Hence 



2 Et 



11 R{i 'j )d " O 2 . 

J<i</<M 



(Ei<y<i rf;i+2dii+Ei</<M */) 



. i=l 



(3.4) 



<9z" !l ■ ■ • dz" lM 



\l/2 



^ R(i,j)z,Zj 

{l<i,j<M 



('"l! ' ' ' '"AiOl J 



;=0 



I 



[j R(i, /)•'' • (3.5) 



Since we have from the definition of e(m\, . . ., ium) that 



1<!</<M 1=1 

for jd,,} e e(mi, . . . , m M ), (3jJ is derived from d3.2t and d3.5t . □ 

We need a simple representation of the coefficient of each v h i> !2 ■ ■ ■ v if in E [exp(Zj ) • • • exp(Z M )] 
where ■ ■ ■ ,k) e (0, 1, • ■ • , d\ e and Zj, . . . Z M are _£]R((A))-valued random variables 
constructed with Gaussian random variables satisfying tl.lOt . 

For C,M e Z >0 , let K t {M) = {k = fa, ...,k M )e (Z> ) M + • • • + k M = {}■ For 
w = v k ■ ■ ■ v i( e A\ let N w : {0, 1, . . . , d\ x {1, . . . , M} x %>(M) Z> be a function 
such that 

N w (i, j, k) = card ({r | i r = i for k x + • • • + + l<r<fci + -- - + kjYj . 

Theorem 3.2 Let w = i?;, t>, 2 • ■ ■ c,, e A* and n w (i) = card ({;' e jl, . . . , €] \ i, = i\) for 
i — \,...,d. Then the coefficient ofw, C{w), in £[exp(Zi) ■ ■ • exp(ZM)] becomes as follows: 



Ifn u '(i) is odd for some is {1, . . . , d\, then 

C(w) = 

Ifn w (i) is even for every ie [1, . . . ,d], then 

1 M 

C ^= E £x^rll( c /) 



(3.6) 



k\ \ ■ ■ - kM'- 
k=(h,...,k u )<=<K c (M) ;=i 



P =i 



I 



e(N u '(p l l,k) r - l N w (p r MM)) 



lll<;<;<M^"'/-J l<i<]<M 



(3.7) 
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where Cj and _R !; are real numbers defined in dl.lO|l . 

Proof In the case where n w (i) is odd for some i e {1, . . . , d\, < l3.6t is directly derived 
from fTTQl 

We therefore consider the other case. By the Taylor expansion of exp (Zi) • • • exp (Z M ), 
we have 

f d \ h ( d 

c\v o + ^ S\Vi ••■ c M v + ^S l M Vi 



°° 1 

£[exp(Zi)---exp(Z M )] = Y r-. t—.E 

h,...,k M =0 

Hence 

C(w) = (E [exp (Zi) ■ ■ ■ exp (Z M )] , w) 



i=i ; 



E 



k=(/d MeftHM) 



I 



k=(ti ImW(M) 



h\---k M \ 



-.11 q iA 'l C^l +1 q^+fcj c t *l- t ~ ,+ *M-l +1 C% + '' + *iV 

V" 5 1 5 2 " ' S 2 " ' S M M 



*A1 



( Cl r (o ' u) • • • (cm)^™ {s\f (14 ' k) ■ ■ • (s^f aM ' k) 



(3.8) 



From the definition of S' „ 



C(w)= L fcl ....fc M . ll( c /) ll E ( s i) •••( s m) 

k=(/ Cl k M )<?Ke(M) l ' M ' ;=1 p=l L 

Applying j3.lt from Lemma l3Tl we obtain <3.7| |. 

On the other hand, the value of the coefficient of each • • ■ c, f in 
j„, (exp (v + (1/2) Yh=\ tf)) can be obtained by the following proposition. 

Proposition 3.3 Let A = [vq, V\V\, V2V2, ■ ■ ■ , VdVd] C A*. Then 

d \ 



exp 



a>i,...,i» ( £/t° 



(3.9) 



Therefore, taking {Sy} . to equate J3.6I I or d3.7t with d3.9l l for w = v h u,- 2 • • ■ c„ 

with \\w\\ < m, we can construct Z\, ... , Zm- 

For m — 5, we take M = 2 to obtain solvable simultaneous equations which in 
fact become the following five: 

1 1 
ci + c 2 = 1, -(ca&n + C2-R22) + #12 = 2' 



1 11 
-(ciR u + c 2 K 22 ) + -ci(i?i2 + R22) = -, 

2 4 

1 11 
-(ciR n + C2R22) + tC 2 (J?h + #22) = 7, 
2 4 

^(R z n + R\ 2 ) + iR u (R« + -R22) + ^11^22 = g. 



(3.10) 
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The solution is <1.15l . Since we let [S'.}j = i r .. ) i i j = i r .ja be the Gaussian system, such 
random variables can be constructed. 

Remark 3.4 If we let m = 5, then M must be at least two. 
4 The Runge-Kutta method 

We begin by briefly introducing the tree theory following 1 3 |, 1 6 \, and |7J. For details 
of the Runge-Kutta method, see GO, Q, and (26l . 

All trees introduced here are called directed or rooted trees in the literature listed 
above. 

Definition 4.1 A labelled tree t is a pair of finite sets ( V(t), £(t)) that satisfies the following 
conditions: 

(1) V(t) c Z, V(t) + 0, and £(t) c {(x, y) e V(t) x V(t) : x < y). 

(2) For each x e V(t), if{x, y) e £(t) and (x' , y) e E(t), then x = x' . 

(3) For two distinct elements x, y e V(t), one of the followings holds: 

(i) Tfere exists a path from x to y. 

(ii) There exists a path from y to x. 

(iii) For some z e ^(t) \ {x, y), there exist paths ztox and z to y. 

Here a path from p\ to pi is a sequence (pi, p-i), (j>2, P3), ■■■ , (pi-i,pi) of elements o/£(t). 

An element ofV(t) is called a vertex oft and that o/E(t) is called an edge of t. 

A particular labelled tree %i is that with card {V(j( )) = 1 and E(t( ) = 0. 

For a labelled tree t = (V(t), E(t)), let r(t) be card (V(t)). We de/ine T as fte set of all 
labelled trees. 

Proposition 4.2 For eacft t = ( V(t), E(t)), t^ere exists a unique vertex r e V(t) such that 
for any x 6 V(t) \ {r\, there is a path from r to x. 

Such a vertex r is called the root of t. Here, Tf consists of only the root. 

Definition 4.3 For i = 1, . . . , n, let t, = ( V (t,) , E (t,)) e T foe suciz f/wf V (t,) n V (ty) = 
ifi * ;'. T/jen [ti • • ■ t„] is defined as t = (V(t), E(t)) e T such that 

v(t) = wuv(ti)u-uy(t„) 

E(t) = {(r, n), . . . , (r, r„)} U E (t x ) U • • • U E (t„) 
where each r{ denotes the root of t, and r = minjri, . . . , r n } — 1. 
Remark 4.4 For ti, . . . , t„ e T, we have that 

[tl • • • t„] = [ta)(i) • • • ta>(n)\ 

for any permutation a) € S„. 

Definition 4.5 Let t, = {V (t,) , E (t,)) e T/or i = 1, 2. We say fZwf ti and t2 are isomorphic, 
written as ti~t2, if there exists a bijection co : V(t\) — > V fa) such that (x,y) e E(ti) if 
and only if(u>(x),u>(y)) € E (t2). 

In particular, when ti~t2 and V (ti) = V (t2), that is, a> is a permutation, we say that ti 
and \.2 are equivalent and write ti ~ t2- 
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Proposition 4.6 Both ~ and ~ are equivalence relations. 

Proposition 4.7 Let t, = (V (t,) , £ (t,)) e T and u, = {V (u,) , E (u,)) e J for i = 1, . . ., n. 

Suppose that t,~u,/or i = 1, . . . , n and that 

V ft) n V (t y ) = and V (uj) n V (u ; ) = 

if z # ;'. T«en 

[ti---t„]~[ui---u„]. 

Definition 4.8 We de/me T = T/~. An element t e T is called a non-labelled tree. For a 
labelled tree t e T, |t| denotes the corresponding non-labelled tree t e T. 

Then, from Proposition l4.7[ the following result can be derived. 
Proposition 4.9 Under the same condition as Proposition \4.7\ 

|[ti---t„]| = |[U!---U H ]| 

holds. 

By virtue of Proposition 14.91 we can define a non-labelled tree t = [t\ ■ ■ ■ t n ] for 
t\,...,t„ € T as |[ti • • • t„]| where t, e T is a representative labelled tree such that 
| tj| = tj. In particular, we let t = \it \. 

Proposition 4.10 For any t e T \ (t), there exist t\,...,t n € T such that t — [ti ■ ■ ■ t n ]. 

Moreover 

[h ■ ■ ■ t„] = [f ffl (i) • ■ ■ 

for any permutation a) € S„. 

Here, [t™ 1 ■ ■ ■ C'] denotes [h ■ ■ ■ h ■ ■ ■ t n ■ ■ ■ t„] where t, 6 T for i = 1, . . . ,n. 

Definition 4.11 (1) For t = (V(t), E(t)) e T, we define a:T — > Z>i, r:T — > Z>i, and 
a : T — > Z>i by 

a{t) = card ({u e T | u ~ t where t e T is a representative element such that |t| = tfj 
r(t) = card (V(t)) 

(l ift = t 

[TlLmioWF ift = \e?-~t?],l>\ 

where A = {aij). ._ ■ We notice that a is well-defined because a denotes the number 

of ways a tree may be labelled. 
(2) Let be the set of K x K real matrices. We inductively define derivative weights 
Q:TxJ{ — >B.fori = l,...,Kby 



a(t) 



In addition, we define the elementary differentials D : C~(R N ;]R N )xr — > C"(R N ;]R N ) 
as follows: 



D(W, t)(x) 



W(x) ift = t, 

W®(x) (D(W, fi)(x), D(W, f 2 )(x), . . . , D(W, f ; )(x)) ift = [ht z --- 1,], I > 1. 

(4.1) 
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Let y( W, s) be a solution to an ODE 



-y(W /S ) = W(y(W,s)), y(W,0) = y 



(4.2) 



where W e C ; "J° (R N , R N ) and yo 6 R N . Then we have the following lemmas essentially 
proved in 0, pp. 139-145. 

Lemma 4.12 For m e Z>\, 

y (m) (W,s)= ^ D(W,|t|)(y(W,s)). (4.3) 

teT„, 

Let T,„ = e T : r{t) = m\ an d T< m = \Jn=o T >« for zrz > with T = 0. 

Let A fx e 3i denote A in 41.161 for the explicit Runge-Kutta method. Then 
Yi(W, s) is definitely determined by A ex with «« = if z < / and so Y(yo; W, s) can be 
constructed with b and Y,(W,s) as both seen in 41.16b . 

Lemma 4.13 Let m > 1. If there exists a constant C m > swe/z £to 



Y,(W,s) 



y° + E * r(i) §§*W)(yo) 



a(f) 



/or z = 1,...,K, then there exists a constant C m+ \ 



l=l,...,m-\ 
t=[h~ti\aT< m 



< C m s m \\W\\™, 



< C m+1 s m+1 ||W||^. 



(4.4) 



(4.5) 



Applying these lemmas to evaluations of the solution to 44.2ft and the Runge- 
Kutta method fTTl6b , we obtain the following result. 

Theorem 4.14 For y satisfying 44. 2b , there exists a constant C m+ \ 



exp (sW)(y ) 



yo- 



+ .E ^«(WV,0(y ) 



< C m+ is 



ni-r-li|TA7l|»'+l 



(4.6) 



On the other hand, for the Runge-Kutta method 41.161 1 there exists a constant C', such 
that 



Y(y ;W,s)- 



,r(t) 



I 



y° + E iW) E fci 1 II C; A) D( n t} ^ 

/=!,.. .,m-l ^ ' 1=1 k=l 



f=[t r -t,]6T<„ 



< c; +1 s m+1 iiwn^ 



We say that (A, b) satisfies wz-th-order conditions if 



(4.7) 



(4.8) 



r(t)\ a(t) 

for alH = [ti---fj]er< ffl . 

From Theorem l4.141 the following result can be directly derived. 

Theorem 4.15 Suppose that (A, b) satisfies the mth-order conditions 44.81 1. Let g(W)(yo) = 
Y(yo; W, 1) where Y(yo; W, 1) is the Runge-Kutta method defined in 41.16b . Then 

g e IS(m). (4.9) 
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5 The new simulation scheme and Corollary 11.41 

Corollary 11.41 indicates the new implementation method of the new higher-order 
scheme proposed by Kusuoka in 1 15 [ and 1 16 [. 

This implementation method seems to be distinct mainly because it has two 
advantages. One is that the approximation operator can be obtained by numerical 
calculations if the Runge-Kutta method is applied to the calculation of each exp (Z ; ) 
whereas the tediousness in symbolical calculations of the operator might be an 
obstacle for practical application, which can be observed in 1 18 |, |24J, and 1 28 1 . The 
other advantage is that the partial sampling problem discussed in [18 1 and [24J can 
be resolved by using quasi-Monte Carlo methods. More precisely, the following two 
points make an effective use of the Low-Discrepancy sequences, which are essential 
to quasi-Monte Carlo methods (| 22) ): 

- In this implementation, S'. can be taken to be a continuous random variable. 

- The scheme itself is characterized by the need for a much less number discretiza- 
tion time steps, which leads to a reduction in the number of dimensions of the 
numerical integration. 

6 Application 

In this section we present a numerical example in order to illustrate the imple- 
mentation method proposed in Corollary 11.41 and compare it with some existing 
schemes. 

6.1 Simulation 

Let X(t,x) be a diffusion process defined by i ll. It . The most popular scheme of 
first order is the Euler-Maruyama scheme, which is shown in [14) and |32"1 . for an 
arbitrary C 4 function / 

||E [/ (X»!)] - £ [/ (X(l, x))]\\ < C/l (6.1) 

where X' EM ''"i denotes the Euler-Maruyama scheme approximating X(t, x). We note 
that this inequality holds for measurable / if { V-, ) (=1 d satisfies some more conditions 
(□)■ 

The construction of a higher-order scheme is based on the higher order stochas- 
tic Taylor formula (|8||14|). When the vector fields {V,}^ =0 commute, higher-order 
schemes can be simplified to a direct product of one-dimensional problem as seen 
in (141 . In contrast, for non-commutative {Vi}f_ 0/ the acquisition of all iterated in- 
tegrals of Brownian motion is required, which is very demanding. This is done 
m Q5 1,| 20 1,| 30 1, 1 31 1 and |18| and generalized as the cubature method on Wiener 
space (Ell). 

Once a pth-order scheme {X( ord P'' M fc/„}fc = o 7 ...,,, is obtained and expanded with some 
constant Kf as 



E[/(xW a )] - E[/(X(l,x))] = JC 7 1 + O(^), 



(6.2) 
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the (p + l)th-order scheme can be derived as 

_3L_ E [/(xtordp),^)] _ _L_ E [/(xN^)] . (6.3) 

This boosting method is called Romberg extrapolation and is shown to be applicable 
to the Euler-Maruyama scheme under certain conditions (|32|). 

The simulation approach must be followed by the numerical calculation of 
. However, when n x d is large, it is practically impossible to proceed 
with the integration by using the trapezoidal formula and so we fall back on the 
Monte Carlo or quasi-Monte Carlo method (|22 1). Here we make only a few remarks 
on each method. For a more detailed analysis, see 1251. 

Remark 6.1 As long as we use the Monte Carlo method for numerical approximation of 
E[f (X(l, x))], the number of sample points needed to attain a given accuracy is independent 
of the number of the dimensions of integration, namely both the number n of partitions and 
the order p of the approximation scheme. 

Remark 6.2 In contrast to the Monte Carlo case, the number of sample points needed for the 
quasi-Monte Carlo method for numerical approximation of E[f(X(l,x))] heavily depends 
on the number of the dimensions of integration. The fewer the dimensions, the fewer the 
samples that are needed. 



6.2 The algorithm and competitors 
6.2.2 The algorithm of the new method 

We take the algorithm which is proposed in Theorem 11.61 and Corollary 11.41 with 
u — 3/4. From Corollary 11.41 we can implement the second-order algorithm with 
a numerical approximation of exp (Z,) of at least fifth-order Runge-Kutta method 
because the order m for an integration scheme attained by Z\ and Z 2 is five and so 
the order of the new implementation method becomes two. As a result of the same 
argument it can be shown that at least seventh-order explicit Runge-Kutta method 
has to be applied to the approximation of exp (Z,) when we boost the new method to 
the third order by Romberg extrapolation. Details of these Runge-Kutta algorithms 
used here are given in the Appendix. 

6.2.2 Competitive schemes 

There there are numerous studies on the acceleration of Monte Carlo methods f fTTI '). 
We choose for the following reasons only the crude Euler-Maruyama scheme and 
the algorithm introduced in 1 25 [, which we will refer to in the remainder of this paper 
as N-V method, both with and without Romberg extrapolation, as competitors: 

(i) Only these two schemes can be recognized as being comparable to the new 
method, since they are model-independent. 

(ii) Almost all variance reduction techniques and dimension reduction techniques 
applicable to the Euler-Maruyama scheme are also applicable to the new 
method. 
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6.3 Numerical results 

We provide an example on financial option pricing in the following part of this 
paper. 

6.3.2 Asian option under the Heston model 

We consider an Asian call option written on an asset whose price process follows 
the Heston stochastic volatility model. Comparison with the N-V method will also 
be given as well from the result shown in |25 [. 

The non-commutativity of this example should be noted here. 
Let Y\ be the price process of an asset following the Heston model: 



where x = (x.\,x-i) e (R>o) 2 , (B 1 (f),B 2 (f)) is a two-dimensional standard Brownian 
motion, -1 < p < 1, and a, 6, \i are some positive coefficients such that 2ad - f} 2 > 
to ensure the existence and uniqueness of a solution to the stochastic differential 
equation (| 9 1). Then the payoff of Asian call option on this asset with maturity T and 
strike K is max (Ys(T, x)/T - K, 0) where 



Hence, the price of this option becomes D x E [max {Y^(T, x)/T - K, 0)] where D is 
an appropriate discount factor that we do not focus on here. We set T — 1, K — 1.05, 
H = 0.05, a = 2.0, p = 0.1, 6 = 0.09, p = 0, and (x 1 ,x 2 ) = (1.0,0.09) and take 



that is obtained by the new method with Romberg extrapolation and the quasi- 
Monte Carlo with n = 96 + 48, and M = 8 x 10 8 where M denotes the number of 
sample points. 

Let Y(t, x) = (Y\(t, x), Yi{t,x), Yj,(t,x)). Transformation of the stochastic differ- 
ential equations l !6.4t and 16.5b gives the following Stratonovich-form stochastic 
differential equations: 




(6.4) 




(6.5) 



E [max(Y 3 (r,x)/T - K, 0)] = 6.0473534496 x 10 -2 




2 



(6.6) 
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Table 1 # of dimensions involved in each method. 



Method 


Number of dimensions 




Euler-Maruyama 


dn 


N-V 


n + dn (n-Bernoulli and (d x n)-Gaussian) 


New Method 


2dn 



where 

( f (yi, yi, y 3 )) = '(yi - f - f ) , «(0 - y 2 ) - yi) 

Vi ( f (yi, y 2 , y 3 )) = \yi VyI, p£ VyE o) ( 6 - 7 ) 
^2 ( f (yi, yi, yi)) = '(o, £ ^{l- P *)y z , o) . 

6.3.2 Dimensions of integrations 

As mentioned in Remarks l6. 1 | and |6T2l the dimensions of integrations in these meth- 
ods affect the quasi-Monte Carlo method. The relation among d: the number of fac- 
tors, n: the number of partitions, and the dimensions of integration of each method 
can be summarized as in Table [U 



6.3.3 Discretization Error 

The relation between discretization error and the number of partitions of each 
algorithm is plotted in Figure [T] We can observe from this figure that for 10" 4 
accuracy the new method with Romberg extrapolation takes the minimum number 
of partitions as n — 1 + 2 whereas n — 16 for the Euler-Maruyama scheme with the 
extrapolation. Even without the extrapolation, the new method attains that accuracy 
with n — 10 while the Euler-Maruyama scheme takes n — 2000. Moreover, it may be 
said that the N-V method shows slightly worse performance than the new method. 

6.3.4 Integration Error 

Looking at Figure EJ we can compare convergence errors of respective methods 
for each number of sample points, M. For the Monte Carlo case, 2a of 10 batches 
is taken as convergence error while for the quasi-Monte Carlo method, absolute 
difference from the value to be convergent is considered. For 10~ 4 accuracy with 
95% confidence level (2a), M = 10 8 is taken for the Monte Carlo method. On the 
other hand, if we apply instead the quasi-Monte Carlo method, the new method 
and the N-V method require M = 2 x 10 5 sample points, though M - 5x 10 6 has to 
be taken for the Euler-Maruyama scheme. 
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Discretization Error and Num. of Partitions 



1e-05 



1e-06 



- \ 


Euler-Maruyama - 
Euler-Maruyama + Romberg - 

Ninomiya-Victoir 
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\ 5e-05 
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'■ ■ 
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1 10 100 1000 10000 100000 

Num. of Partitions 

Fig. 1 Error coming from the discretization 




1000 10000 100000 
Num. of Sample points 



1e+06 



1e+07 



Fig. 2 Convergence Error from quasi-Monte Carlo and Monte Carlo 
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Table 2 #Partitions, #Samples, Dimension, and CPU time required for an accuracy of 10 4 . 



Method 


#Part. 


Dim. 


#Samples 


CPU time (sec) 




E-M + MC 


2000 


4000 


10 s 


1.72 x 10 5 


E-M + Romb. + QMC 


16 + 8 


48 


5x 10 6 


1.27 x 10 2 


N-V + QMC 


16 


32 + 16 


2xl0 5 


4.38 


N-V + Romb. + QMC 


4 + 2 


12 + 6 


2xl0 5 


1.76 


New Method + QMC 


10 


40 


2xl0 5 


3.4 


New Method + Romb. + QMC 


2 + 1 


12 


2xl0 5 


1.2 



6.3.5 Overall performance comparison 

The number of partitions, the number of samples, and the amount of computation 
time required for 10" 4 accuracy for each method are summarized in Table [2] CPU 
used in this experiment is Athlon 64 3800+ by AMD. 

Since the amount of time required to carry out the calculation for each sample 
point is proportional to the number of partitions, the total time spent on calculations 
is proportional both to the number of partitions and to the number of samples. We 
can see from the Table [2] that the speed of the new method is approximately 100 
times faster than that of the Euler-Maruyama scheme when Romberg extrapolation 
and quasi-Monte Carlo are applied to each. Even when the extrapolation is not 
applied, the new method enables calculations some 37 times faster than the Euler- 
Maruyama scheme with Romberg extrapolation and quasi-Monte Carlo method. 
This fact shows that the reduction in the number of partitions sufficiently compen- 
sate for the slowness of one step of the new method at least in the present study. 

Lastly, Remarks l6.1| and |6T2l should be emphasized to reiterate that the advantage 
of the new method is that it is deeply related to the properties of the quasi-Monte 
Carlo method. 



Appendix: The fifth-order and the seventh-order Runge-Kutta algorithms 

We present here the concrete algorithms of the explicit fifth- and seventh-order 
Runge-Kutta methods applied in Subsection 16.21 The fifth-order method is taken 
from [ 6 [ as follows: 



2 


11 




5 


1 


3 


15 


5 


* 31 = 64' 


«32 = 


= 64' 


a 43 - 2' 


fl51 = 64' 


^ = -64' 


3 


9 




_ 5 


6 


12 


8 


a 53 = -, 


* 54= 16' 


«62 


~~ 7' 


«63 = J, 


#64 ="/ 


«65 = J, 






Oij 


= 


otherwise, 







7 . 32 12 32 7 
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The seventh-order method is taken from Q as follows: 

_ 1 _ 1 _ 1 _ 3 148 150 56 

#21 — #32 — 7/ #41 — „/ #43 — „/ #51 — 7777/ #53 — 7777/ #54 



fi , "32 -g, -41 "g, -43 - g , "51-^, ^" 1331 , -54- ^ , 

404 170 4024 10648 2466 1242 



#61 - "243, #63 - "27, #64 - ^ «65 " W > #71 - ^J, #73 - 343 , 

19176 51909 1053 5 96 1815 



a?i ~ "16807' * 75 " "16807' ° 76 " 2401' " 81 " 154' 084 " 539' " 85 " "20384' 
405 49 113 195 32 29403 



a86 ""2464' fl87 "Tl44' fl91 ""l2"' 893 " "a"' 094 " T' " 95 " ^58T' 



#96 = -777/ #97 = 777,77' #98 - 77^ a 'j - otherwise, 



729 _ 1029 _ 21 

"512' 097 " 1408' 098 "16' 

, , ( , 32 1771561 243 16807 77 11 x 
o=000 



105 6289920 1560 74880 1440 70/ 



References 



1. Bally, V., Talay, D.: The law of the Euler scheme for stochastic differential equations I. Con- 
vergence rate of the distribution function. Probability theory and related fields 104, 43-60 
(1996) 

2. Bayer, C, Teichmann, J.: Cubature on Wiener space in infinite dimension. preprint: 
arXiv:8712.3763[math.PR] (2007) 

3. Bollobas, B.: Graph Theory: An introductory Course. Springer Verlag (1979) 

4. Bourbaki, N.: Elements de Mathematique, Croupes et Algebres de Lie, Chapitres 2 et 3. Her- 
mann, Paris (1972) 

5. Burrage, K., Burrage, P.M.: General order conditions for stochastic Runge-Kutta methods for 
both commuting and non-commuting stochastic ordinary differential equation systems. Appl. 
Numer. Math. 28(2-4), 161-177 (1998) 

6. Butcher, J.C.: The Numerical Analysis of Ordinary Differential Equations. John Wiley & Sons, 
Chichester (1987) 

7. Butcher, J.C.: Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, 
Chichester (2003) 

8. Castell, E: Asymptotic expansion of stochastic flows. Probability theory and related fields 96(2), 
225-239 (1993) 

9. Feller, W.: Two singular diffusion problems. Annals of Mathematics 54, 173-182 (1951) 

10. Filipovic, D., Tappe, S., Teichmann, J.: Jump-diffusions in Hilbert spaces: existence, stability and 
numerics, preprint: arXiv:881S. 582 3 [math. PR] (2008) 

11. Glasserman, P.: Monte Carlo Methods in Financial Engineering. Springer verlag, New York 
(2004) 

12. Hairer, E., Lubich, C, Wanner, O: Geometric Numerical Integration: Structure-Preserving Al- 
gorithms for Ordinary Differential Equations, 2nd ed. Springer Verlag (2006) 

13. Ikeda, N, Watanabe, S.: Stochastic differential equations and diffusion processes. North Hol- 
land/Kodansha (1981) 

14. Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic Differential Equations. Springer 
Verlag, Berlin (1999) 

15. Kusuoka, S.: Approximation of Expectation of Diffusion Process and Mathematical Finance. 
In: T. Sunada (ed.) Advanced Studies in Pure Mathematics, Proceedings of Final Taniguchi 
Symposium, Nara 1998, vol. 31, pp. 147-165 (2001) 

16. Kusuoka, S.: Approximation of Expectation of Diffusion Processes based on Lie Algebra and 
Malliavin Calculs. Advances in Mathematical Economics 6, 69-83 (2004) 

17. Kusuoka, S.: Kusuoka Scheme and Gaussian type approximation. Presentation at "Mathemat- 
ical Finance Seminar in Graduate School of Mathematical Sciences The University of Tokyo 
(l/June/2005)" (2005) 

18. Kusuoka, S., Ninomiya, S.: A new simulation method of diffusion processes applied to Finance. 
In: J. Akahori, S. Ogawa, S. Watanabe (eds.) Stochastic processes and application to mathe- 
matical finance, Proceedings of the Ritsumeikan International Symposium, pp. 233-253. World 
Scientific, Singapore (2004) 



25 



19. Lapeyre, B., Pardoux, E., Sends, R.: Methodes de Monte-Carlo pour les equations de transport 
et de diffusion (Mathematics and Applications 29). Springer Verlag, Berlin (1998) 

20. Liu, X.Q., Li, C.W.: Weak approximation and extrapolations of stochastic differential equations 
with jumps. SIAM Journal on Numerical Analysis 37, 1747-1767 (2000) 

21. Lyons, T., Victoir, N.: Cubature on Wiener Space. Proceedings of the Royal Society of London. 
Series A. Mathematical and Physical Sciences 460, 169-198 (2004) 

22. Niederreiter, H.: Random Number Generation and Quasi-Monte Carlo Methods. SIAM (1992) 

23. Ninomiya, S.: A new simulation scheme of diffusion processes: Application of the Kusuoka 
approximation to Finance Problems. Mathematics and Computers in Simulation 62/3-6, 479- 
486 (2003) 

24. Ninomiya, S.: A partial sampling method applied to the Kusuoka approximation. Monte Carlo 
Methods and Applications 9, 27-38 (2003) 

25. Ninomiya, S., Victoir, N.: Weak Approximation of Stochastic Differential Equations and Appli- 
cation to Derivative Pricing. Applied Mathematical Finance 15, 107-121 (2008) 

26. Rossler, A.: Runge-Kutta Methods for the Numerical Solution of Stochastic Differential Equa- 
tions. Shaker Verlag GmbH (2003) 

27. Rumelin, W.: Numerical treatment of stochastic differential equations. SIAM Journal on Nu- 
merical Analysis 19(3), 604-613 (1982) 

28. Shimizu, M.: Application of the Kusuoka approximation with Tree Based Branching Algorithm 
to pricing interest-rate derivatives with the HJM model. Master thesis, Imperial College of 
Scienece, Technology, and Medicine (2002) 

29. Strichartz, R.E.: The Campbell-Baker-Hausdorff-Dynkin Formula and Solutions of Differential 
Equations. Journal of Functional Analysis 72, 320-345 (1987) 

30. Talay, D.: Second-order discretization schemes of stochastic differential systems for the compu- 
tation of the invariant law. Stochastics and Stochastics Reports 29, 13-36 (1990) 

31. Talay, D.: Simulation of Stochastic Differential Systems. In: P. Kree, W. Wedig (eds.) Probabilistic 
Methods in Applied Physics, LNP 451, pp. 54-96. Springer Verlag (1995) 

32. Talay, D., Tubaro, L.: Expansion of the global error for numerical schemes solving Stochastic 
Differential Equations. Stochastic Analysis and Applications 8, 483-509 (1990) 



